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Quantum effects of plasmonic phenomena have been explored through ab-initio studies, but only for exceed¬ 
ingly small metallic nanostructures, leaving most experimentally relevant structures too large to handle. We 
propose instead an effective description with the computationally appealing features of classical electrodynamics, 
while quantum properties are described accurately through an infinitely thin layer of dipoles oriented normally to 
the metal surface. The nonlocal polarizability of the dipole layer is mapped from the free-electron distribution 
near the metal surface as obtained with ID quantum calculations, such as time-dependent density-functional 
theory (TDDFT), and is determined once and for all. The model can be applied to any system size that is tractable 
within classical electrodynamics, while capturing quantum plasmonic aspects of nonlocal response and a finite 
work function with TDDFT-level accuracy. Applying the theory to dimers we find quantum-corrections to the 
hybridization even in mesoscopic dimers as long as the gap is sub-nanometric itself. 


Nanoscale metallic structures receive considerable atten¬ 
tion due to plasmon phenomena beyond classical electrody¬ 
namics [1-5]. In deep nanoscale structures, quantum effects 
will manifest themselves [6]. Size-dependent plasmon ener¬ 
gies of nanoscale metallic particles, an accustomed topic in 
cluster science and surface physics [7-11], has been revisited 
with new insights [1, 2]. Dimers hosting Angstrom-scale 
gaps are being investigated intensively in search of quantum¬ 
tunneling effects [3,4], though interpretations highlight both 
quantum tunneling [12] and semi-classical dynamics [13]. 
Time-dependent density-functional theory (TDDFT) allows 
explorations into the quantum regime [14, 16, 17], albeit 
limited to atomic-scale systems. 

Here we focus instead on the mesoscopic scale that is 
typical for state-of-the-art experiments: too small to obey 
the classical local-response approximation (LRA), but too 
large for an efficient quantum treatment. Plasmonics at 
this scale has so far been described with some success by 
semi-classical hydrodynamic Drude models (HDM) that ex¬ 
hibit nonlocal response [5, 13, 18-21] and by the quantum- 
corrected model (QCM) for gaps between two metals [12]. 
It is however desirable to look beyond these models. 

In this Letter, we present the projected-dipole model 
(PDM), a theory that describes optical properties of meso¬ 
scopic plasmonic systems with TDDFT level accuracy. Their 
quantum effects are captured by a zero-thickness projected- 
dipole layer (PDL) onto the surface of the metal, see Fig. 1. 
The idea is partly inspired by Ref. 19, where a finite- 
thickness local layer represents nonlocal response within 
the HDM. Here, we go beyond the semi-classical HDM 
and reproduce quantum phenomena of a microscopic theory 
of choice, TDDFT. A crucial computational advantage of 
the PDL is that its inclusion can be elegantly absorbed into 
boundary conditions at the interface. 

Projected-dipole layer .— We first introduce the PDL, 
the key theoretical concept of our approach. To repro- 



FIG. 1 (color online). In the projected-dipole model (PDM), light- 
driven quantum-plasmonic response in a metallic nanostructure 
(left) is represented in classical local-response electrodynamics 
by an infinitely thin layer of dipoles that point perpendicularly 
to the surface (right). The loupe (left) highlights the spill-out of 
the microscopic electron distribution inducing a dipole moment 
proportional to the Feibelman parameter de¬ 


duce a quantum calculation of the reflectivity of a planar 
metal-air interface at x = 0 within a classical model, we 
add an infinitely thin layer of dipoles to the corresponding 
classical interface. The electric response of this PDL is 
PpDL = 6otrnn • Eair(x = 0), where n is the normal vector 
pointing outward the metal surface and Eak is the electric 
field in the air region. The polarization field Ppdl is a projec¬ 
tion of the electric dipole onto the surface normal, with a yet 
unknown polarizability a that we identify uniquely below, 
after further motivating our model. The semi-infinite metal 
at X < 0 that is in the electrostatic limit driven by an external 
electric potential 0ext = where s = yy + zz, exhibits 

an induced electron density p(x) that within the LRA is a 
Dirac delta function qiS(x). The induced electric potential 
outside the metal equals I2k\\. By contrast, p{x) 

spreads across the geometrical boundary in quantum theo¬ 
ries. The corresponding induced electric potential outside 
the metal consequently has a multipole expansion that starts 
as (-^^2e“^"^/2/:||)[l -F k\\dc + Oik^dl)]. Here the first and sec¬ 
ond terms in the bracket represent the monopole and dipole 
contributions, respectively, with q 2 = f dxp(x) the total in¬ 
duced electron charge and dc = f dx xp(x)lq 2 the Feibelman 
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parameter, i.e., the centroid of the induced charge. Clearly, 
the leading-order quantum effect is to induce a dipole mo¬ 
ment proportional to dc in the surface region as indicated in 
Fig. 1. This motivates our choice of a PDL within classical 
electrodynamics to model quantum effects. 

Planar-surface polarizability .— Next we outline in more 
detail how to identify the polarizability a of the PDL. Within 
TDDFT, we solve p(v, cS) self-consistently from 


p(v, (jj)=J dx' [S(v, v', oj)p(x, oj)-exKs(x, x\ u;)0ext(-^', oj)]. 

Here, Xks represents the Kohn-Sham response function and 
S accounts for Coulomb and exchange-correlation interac¬ 
tions [22]. For simplicity we employ the jellium approxi¬ 
mation of the positive-ion background and the local-density 
approximation of the exchange-correlation potential. 

We calculate the TDDFT reflection amplitude /?ks(^ii, oj) = 
fdxp(x, 6G))e^ii^, and require that the classical PDL prob¬ 
lem gives the same value, which gives 


a(cL),kii) 




( 1 ) 


where is the bulk permittivity of the metal. To interpret 
this result, we first consider its long-wavelength (or k,, ^ 0) 
limit OfL. Since P^s in this limit equals - 1][1 

k|,r/,(a;)]/[eLRA(a;)+1 -(eLRA(^)_ [23], we obtain 

which indeed shows 
a direct relation of the PDL polarizability with dc as our 
intuitive expectation. The latter parameterizes two important 
quantum effects: electron spill-out via Re[Jc] and electron- 
hole excitations via Im[Jc] [7]. Within the HDM, dc becomes 
-l/^L (“spill-in”, proper for noble metals [2, 5] but not 
for simple metals [16, 17]) in terms of the Thomas-Fermi 
wavenumber [18], and becomes (6^^^ - l)/(6j(;^^kL), 
which has the same effect as the local dielectric layer in 
Ref. 19. The standard HDM cannot describe electron spill- 
out or electron-hole excitation, though; only recent HDM 
extensions do include spill-out [20, 24, 25]. 

The plasmonic response of Sodium exhibits an interplay 
between nonlocal response and quantum spill-out due to a 
finite work function, and illustrates the potential of our PDM. 
Fig. 2(a) shows reflection spectra off a planar Sodium-air 
interface for different a standard TDDFT calculation. The 
peak around 4 eV corresponds to the surface plasmon (SP) 
resonance, which first shifts to the red and then to the blue as 
k|| increases, while the width shows a noticeable broadening 
due to enhanced electron-hole excitations. Fig. 2(b) depicts 
OfL, with around 4.8 eV a sign change of Re[QfL] and a corre¬ 
sponding peak in Im[a'L]. This multipole surface plasmon 
(MSP) resonance [8], with dipole-like charge distribution at 
the surface, is barely visible in Fig. 2(a) due to the nearby 
stronger SP (monopole) resonance. Fig. 2(c) shows an in¬ 
creased /:||-dependence of a for larger energies, so is a 
good approximation of a only at low energies. 



FIG. 2 (color online), (a) TDDFT reflection spectra of a Sodium 
surface for wave incidence with different k\\ (in units of ^p)- (b) 
Long-wavelength polarizability ul. (c) Magnitude of polarizability 
a normalized by Up. 


Arbitrary structures. — With the flat-layer polarizability 
k||) defined in Eq. (l)-the input parameters of the theory, 
we here generalize it to an arbitrary metal structure with a 
locally curved surface. Thus, we seek the polarizability 
in real space, which for the planar surface is given by the 
inverse Fourier transform of which assumes the 

nonlocal form apico, |s-s'|) that only depends on the distance 
between the points on the surface. As observed in Fig. 2(c), 
a is nearly constant for k,, < kc with kc ~ O.lkp. Thus the 
nonlocal range of the flat-layer is of the order l/kc ^ 1 nm. 
For a general curved surface as long as the local curvature 
does not exceed 1 nm“\ we make the approximation that 
the polarizability of a curved surface has the expression of 
a planar surface, with the distance on the curved surface 
defined as the length of the geodesic [26] . We thus obtain 

a(a), s, s') = Xj Is - s'lg). (2) 

The sum over g runs over all geodesics between s and s', 
with length |s - s'|^. When g = I and the geodesic is a 
straight line, Eq. (2) simplifies to the planar-surface case. 

By using Eq. (2) for a, the response of an arbitrary struc¬ 
ture including quantum effects can be predicted within classi¬ 
cal electrodynamics. With retardation included, the presence 
of the PDL is equivalent to the boundary condition 

R||, air(s) “ E||,metal(s) — ds (XicO, S, S )n(s ) * Eair(s ), 

( 3 ) 

where the subscript ’’H” represents the parallel components 
of the electric fields [27] . This boundary condition can be 
conveniently adopted in numerical computations. If large-k,, 
fields do not contribute to the dynamics significantly or if one 
is just interested in the leading-order quantum effects, a can 
be approximated by the local form QfL(u;)h(s - s'), as can be 
derived by replacing a{k\\, co) with a^ico). The boundary con¬ 
dition simplifies to E,,, air - E,,, metal = (n • Eair). Further, 

beyond the local approximation, with the observation that 
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the nonlocal polarizability is a short-range function, we can 
make the shortest-geodesic approximation, i.e., just keeping 
the term of the shortest geodesic length in Eq. (2) [27]. 

In the following we focus on plasmonic structures with 
ID translational invariance, in the z-direction say, which in¬ 
clude extended ID gratings and arbitrarily shaped nano wires. 
This simplifies the task of finding geodesics. After Fourier 
transformation in the z-direction with wavenumber the 
polarizability of the PDL in Eq. (2) becomes 
+ 00 

a((jtj,kzJ,r) = ^ a(kn, ~ ^ ^. 

n=-oo 

Here L is the total arc length of the boundary in the v, y- 
plane, I g [0, L] is the local coordinate of the boundary 
with Cartesian metric, and kn = [(InnlLf + In this 

picture, fields propagate along the curved boundary as if it 
were a planar interface with artificial periodicity L, so that a 
is spanned by the associated Fourier components. 

Benchmarks against TDDFT — We will confirm the valid¬ 
ity of the PDM for a Sodium cylinder of radius 4.9 nm (i.e. 
curvature assumption behind PDM is fulfilled), excited by a 
plane wave incident normally to the wire axis and with an 
in-plane polarized electric field. The local bulk permittivity 
is described by the Drude model, = 1 - oj^lojioj iy), 
where oj^ = 5.89 eV. For the damping we use y = 0.16 eV, 
and 0.085 eV for the FRA and PDM, respectively, to match 
the resonance width. The extinction spectra predicted by 
FRA, PDM, and TDDFT are presented in Fig. 3(a). The 
FRA and PDM spectra are calculated with the Green’s func¬ 
tion surface-integral method [28], the TDDFT results are 
reproduced from Ref. 16. Since retardation is neglected 
in the TDDFT, the FRA and PDM are also calculated in 
the electrostatic limit. The PDM and TDDFT spectra in 
Fig. 3(a) agree extremely well, and both predict that the SP 
peak around 4 eV is redshifted with respect to the FRA peak. 
Next, we consider a dimer of such cylindrical nanowires, 
separated by a gap of only 0.74 nm. Fig. 3(b) shows a pro¬ 
nounced discrepancy between the FRA extinction spectrum 
on the one hand and the PDM and TDDFT spectra on the 
other, with nonclassical broadening in the latter two. More¬ 
over, as one of our main results we find that the PDM and 
TDDFT spectra agree well, quite surprisingly given the fact 
that the PDM models both metal-air interfaces in the gap 
region in an independent-surface approximation (ISA). 

Quantum effects are known to be negligible for large- 
scale structures, e.g. nano wires above 10 nm radius. The 
multi-scale dimer system with a sub-nanometer gap distance 
introduces an additional scale for confining plasmon fields 
and enhances the quantum effects [14, 15]. Such a setup is 
usually far too challenging for ab-initio studies. However, 
the PDM can be applied to explore the quantum effects of 
such a system. Fig. 3(c) depicts the color map of the extinc¬ 
tion cross section (normalized by wire radius) predicted by 
the PDM for the dimer structure as in Fig. 3(b) with radius 
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FIG. 3 (color online). Extinction properties of Sodium nanowires 
within FRA, PDM, and TDDFT. (a) Cylindrical wire with radius 
4.9 nm and (b) corresponding dimer with a gap of 0.74 nm. (c) 
Dimer strueture with radius varying from 3 nm to 30 nm and 0.74 
nm gap distance, where the extinction cross sections are normalized 
by the nanowire radius, and (d) 20 nm radius case [dash-dotted line 
in (e)]. TDDFT results in (a,b) are reproduced from Ref. 16. 


varying from 3 nm to 30 nm. Retardation effects are now 
also included. To demonstrate the quantum effects of a typi¬ 
cal multi-scale system. Fig. 3(d) plots the extinction spectra 
at 20 nm radius of Fig. 3(c) (dash-dot line), contrasting PDM 
and FRA results. Similarly as in Fig. 3(b), quantum effects 
manifest themselves by shifting and broadening the SP res¬ 
onances. Interestingly, the higher-order modes in this case 
are less affected by quantum effects, e.g., the quadrupole 
mode (around 2.8 eV) and hexapole mode (around 3.2 eV) 
show clear peaks in the PDM spectra, since plasmon fields 
in larger structures are less confined. 

More Gap Effects .— In the PDM, we invoked the ISA and 
thus only included quantum effects related to the single metal 
surface. When the gap approaches the regime set by the dc 
length scale, it also implies that a significant wavefunction 
overlap of the two surfaces should be taken into account. 
An intuitive way to do this is to treat the gap as an effective 
medium, as firstly introduced in the QCM with a local Drude 
permittivity constructed to mimic the dissipation associated 
with quantum tunneling (QT) [12]. To improve our PDM 
method for the smallest gaps, here we present an unambigu¬ 
ous extraction of the effective gap permittivity directly from 
TDDFT calculations, instead of relying on the assumption 
that QT is the cause of all gap-related dissipation. 

Turning from the spectral to the spatial information pro¬ 
vided by TDDFT calculations, we extract the effective local 
gap permittivity 6eff, defined as D/[6oE]. In Fig. 4(a) we 
visualize Eqq for a single interface and for two gap distances, 
of 5A and 0.5A. For the 5A gap case, two separated metallic 
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FIG. 4 (color online). Gap effects for planar Sodium-vacuum- 
Sodium gap structures within TDDFT. (a) Effective permittivity 
for single interface and for 5A and 0.5A gaps, (b) Frequency de¬ 
pendence of the corresponding mid-gap effective permittivities, (c) 
Electron-hole backward scattering proportion (BSP) as a function 
of gap distance and energy. 


systems can be identified as seen from the clear boundaries 
between the negative (metal) and positive (gap) valued re¬ 
gions of Re[6eff], as for the single interface. Furthermore, the 
loss Im[6eff] clearly peaks right at the interfaces where we 
likewise have Re[6eff] ~ 0. In other words, we find that dis¬ 
sipation (associated with electron-hole excitations) occurs 
around the metal surface, in agreement with mesoscopic 
transport theory [29], rather than in the vacuum region itself. 
By contrast, for the 0.5 A gap no definite boundary of the two 
metals can be observed from the equilibrium electron density 
distribution. The induced electron current at the gap center 
is accordingly the classical Drude current. In the optical 
spectra this regime (gaps of 1-2A and below) is character¬ 
ized by a transition from the bonding-dipole plasmon (BDP) 
resonance to the charge-transfer plasmon (CTP). 

Fig. 4(b) shows mid-gap properties, 6gap = 6eff(x = 0), for 
gaps of 5A, 2A, and 0.5A. In general, we anticipate a transi¬ 
tion from non-Drude (quantum) response to Drude (classical) 
response as the gap reduces. Indeed, for the 0.5A gap the 
mid-gap response is well captured by the classical Drude 
model (dashed line) using a plasma frequency determined 
by the mid-gap equilibrium electron density [30]. On the 
other hand, the larger separation of 5A shows a complete 
non-Drude behavior. The 2A gap exhibits almost frequency 
independent Re[6gap], thus constituting the cross-over gap 
size between non-Drude and Drude responses. 

Modeling A-sized gaps .— We will now discuss how in 
the PDM a dimer gap can be treated as an effective medium 
(in spirit close to the QCM [12]), thereby accounting for 
not only quantum effects of both individual metal interfaces 
(as before, employing the ISA) but also for their strong 
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FIG. 5 (color online). Extinction cross sections of a cylindrical 
Sodium nanowire dimer with radius 4.9 nm, for gaps decreasing 
from 5.3 to OA. Comparison of predictions by ERA, TDDFT, PDM, 
and PDGM. TDDFT spectra originate from Ref. 16. 


mutual proximity. The procedure in this Projected-Dipole 
Gap Model (PDGM) is to extract 6gap as before from TDDFT 
calculations of a planar dimer, and to treat the gap as an 
effective medium with permittivity 6gap within the PDM 
framework. To illustrate the accuracy of the PDGM, we 
consider a nanowire dimer of two Sodium cylinders with 
radius 4.9 nm driven by a uniform electric field along the 
gap. In Fig. 5 we illustrate the extinction spectra according 
to ERA, PDM, PDGM, and TDDFT, as the gap distance 
decreases from 5.3 to OA (contact). As expected, not all 
ERA resonances are supported by the more accurate TDDFT 
calculations. The PDM captures the resonance broadening as 
the gap is reduced, with spectra closely resembling those of 
TDDFT. However, the PDM neglects strong gap effects, and 
therefore underestimates the resonance broadening. For the 
same reason the CTP does not appear before physical contact. 
Turning to the PDGM, we find almost perfect agreement with 
the TDDFT spectra (both in terms of resonance positions 
and broadening) even for the smallest gaps. 

Discussion .— The great accuracy of the PDM in Fig. 3 
and of the PDGM in Fig. 5 calls for a discussion of the dissi¬ 
pative quantum processes at interfaces and near gaps. The 
TDDFT calculations account for plasmon damping due to 
electron-hole (e-h) excitations. To interpret our results, we 
distinguish two scattering processes (details in Supplemen¬ 
tary Material [27]): an electron state that propagates towards 
the interface or gap can create e-h pairs in forward- or in 
backward scattering. For a single interface, at energies below 
the work function all dissipation due to e-h creation is associ¬ 
ated with backward scattering, while forward scattering can 
only exist above the work function. By contrast for a dimer 
structure with a narrow gap, forward scattering below the 
work function may also play a role. The PDM assumes all 
gap- and interface-related loss to be due to backward scatter- 
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ing, whereas the QCM calculates all gap-related loss as due 
to QT, i.e. forward scattering [12]. Finally, the PDGM does 
not systematically exclude one of both processes. So what 
dominates e-h creation, forward or backward scattering? 

To make our discussion specific, we consider a vacuum 
gap formed between two parallel planar Sodium surfaces 
driven by a uniform electric field. Fig. 4(c) depicts the e-h 
backward scattering proportion (BSP) [27], defined by the 
ratio between the backward and total scattering rates. The 
BSP exhibits completely different behavior for gap sizes 
above and below 3A. For gaps > 7A, the BSP approaches 
100% below the work function, which explains the validity 
of the PDM in Fig. 3. For smaller gaps down to 3A, the 
BSP continuously decreases. Below the work function (3.06 
eV for Sodium interface) this is due to QT. The latter mech¬ 
anism also explains the abrupt decrease of the BSP as the 
energy exceeds the work function [dashed line in Fig. 4(a)]. 
As the gap shrinks further below 3A, the concepts of QT 
and work function are not well defined, since the vacuum 
gap separating different systems cannot be unambiguously 
defined. Thus, the traditional tunneling picture is not proper. 
This is evident from the observation that in Fig. 4(c) the BSP 
starts increasing when making gaps smaller than 3A, and 
eventually dominates the forward-scattering at low energies. 

In conclusion, we propose a projected-dipole model 
(PDM) for quantum plasmonics, where the complex dynam¬ 
ics of the quantum electron gas is explored with TDDFT 
(or other ab-initio methods) in the vicinity of the surface 
and then subsequently the quantum response is represented 
within the common framework of classical local-response 
electrodynamics as a sheet of dipoles pointing normally to 
the surface. Remarkably, we have shown that this mapping 
captures quantum plasmonic aspects of nonlocal response 
and a finite work function with TDDFT-level accuracy. Our 
work also provides key insight into the possible importance 
of optical tunneling currents in plasmonic dimers. Finally, 
our results have great potential for multi-scale modelling in 
computational plasmonics, where quantum effects can now 
be accurately explored even in large systems, while enjoying 
the efficiency and accuracy of state-of-the-art electromag¬ 
netic solvers for classical local-response electrodynamics. 
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